O paulo vai escrever um monte de ladainha nesta nova versao...

#############################################
#########global variables
uniq2=c(unique(table_patterns[,2]),"total")
table_patterns=read.table("table_patterns.txt",header=TRUE,colClasses="character")
source("functions.R")

# Incluindo um comentario
thres[["thr_0.01"]]=read.table("island3to4_thr_0.01")
thres[["thr_0.02"]]=read.table("island3to4_thr_0.02")
thres[["thr_0.05"]]=read.table("island3to4_thr_0.05")
thres[["thr_0.07"]]=read.table("island3to4_thr_0.07")
thres[["thr_0.1"]]=read.table("island3to4_thr_0.1")
save(thres,file="list of island3to4.RData")


######################################################################################################################################################
#####resampling of combinations of 3 islands from clusters with samples from 4 islands (24.11.2010)
sum_tab_subtr=NULL
aux_problm=NULL
aux=NULL
matrixCl=NULL

compar=cbind(compar,uniq2[1:42])
colnames(compar)="patterns"

for(y in 1:length(thres))
	{
	matrixCl=thres[[y]]
	all_subtr=list()
	all_nw=vector()
	aux=which(names(matrixCl)==names(thres[y]))
	uniq1=unique(matrixCl[,aux])
	patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
	pat1=row.names(patterns)
	pat2=colnames(patterns)
	print(Sys.time())
	for (x in uniq1)   
		{
		subAux=matrixCl[matrixCl[,aux]==x,] 
		NbIsland=length(unique(subAux[,2])) 
		if(NbIsland==4)
			{
			count_3=count_3+1
			m=grep("moor",subAux[,1],value=TRUE)
			t=grep("tahi",subAux[,1],value=TRUE)
			h=grep("huah",subAux[,1],value=TRUE)
			r=grep("raia",subAux[,1],value=TRUE)
			m_out=do.call("comb_3",list(t,h,r))
			h_out=do.call("comb_3",list(m,r,t))
			r_out=do.call("comb_3",list(m,h,t))
			t_out=do.call("comb_3",list(m,h,r))
			all_comb=cbind(m_out,h_out,r_out,t_out)
			for (j in 1:dim(all_comb)[2]) 
				{
				count_tot=count_tot+1
				subtr_nw=NULL # restart subtr_nw
				subtr=subset(my3_phylo4,tips.include=all_comb[,j])
				all_subtr[[paste(x,j,sep="_")]]=subtr
				subtr_phylo=as(subtr,"phylo")
				subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
				subtr_phylo$edge.length=NULL
				subtr_newick=write.tree(subtr_phylo)
				subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
				sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
				if (length(sumary_subtr)==2)
					{
					sumary_subtr=c(sumary_subtr, "pb")
					}
				sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
				patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1		
				} 
			patterns[pat1==x,pat2=="total"]=count_tot
			}
		}
		write.table(patterns,paste("patterns for ",names(thres[y]),"_subset_3of4_",Sys.Date(),".txt"), sep="\t",row.names=TRUE,col.names=TRUE)
	write.table(sum_tab_subtr,paste("summary of table subtree",names(thres[y]),"_subset3of4.txt"))
	save(all_subtr,file=paste("list of all subtrees ",names(thres[y]),"_subset_3of4",Sys.Date(),".RData"))
	print(Sys.time())	
	}	
	
	
###################################
#######version not optimized but working
h_t1=read.table("patterns for  thr_0.01 _subset3of4_ 2010-11-24 .txt",header=TRUE)
colnames(h_t1)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr1=apply(h_t1[,(1:42)],2,sum)
count_rows_thr1=apply(h_t1[,(1:42)],1,sum)
t1=count_cols_thr1
t1=cbind(t1,names(t1))
colnames(t1)=c("thr001","patterns")


h_t2=read.table("patterns for  thr_0.02 _subset3of4_2010-11-26 .txt",header=TRUE)
colnames(h_t2)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr2=apply(h_t2[,(1:42)],2,sum)
count_rows_thr2=apply(h_t2[,(1:42)],1,sum)
t2=count_cols_thr2
t2=cbind(t2,names(t2))
colnames(t2)=c("thr002","patterns")
compar_t1t2=merge(t1,t2)
#compar_t1t2=compar_t1t2[,-2]
#topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
#table_comp=merge(compar_t1t2,topo)

h_t5=read.table("patterns for  thr_0.05 _subset3of4_ 2010-11-25 .txt",header=TRUE)
colnames(h_t5)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr5=apply(h_t5[,(1:42)],2,sum)
count_rows_thr5=apply(h_t5[,(1:42)],1,sum)
t5=count_cols_thr5
t5=cbind(t5,names(t5))
colnames(t5)=c("thr005","patterns")
compar_t1t2t5=merge(compar_t1t2,t5)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-c(3,5)]
table_comp=merge(compar_t1t2t5,topo)
table_comp3of4=table_comp[table_comp$isl3_pat_only=="y",]
write.table(table_comp3of4,"comparison thr001,002 and 005 for 3of4.txt",sep="\t")
	
	
h_t7=read.table("patterns for  thr_0.07 _subset_3of4_ 2010-11-28 .txt",header=TRUE)
colnames(h_t7)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr7=apply(h_t7[,(1:42)],2,sum)
count_rows_thr7=apply(h_t7[,(1:42)],1,sum)
t7=count_cols_thr7
t7=cbind(t7,names(t7))
colnames(t7)=c("thr007","patterns")
compar_t1t2t5=read.table("comparison thr001,002 and 005 for 3of4.txt",header=TRUE)
compar_t1t2t5t7=merge(compar_t1t2t5,t7)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-c(3,5)]
table_comp=merge(compar_t1t2t5t7,topo)
table_comp3of4=table_comp[table_comp$isl3_pat_only=="y",]
write.table(table_comp3of4,"comparison thr001,002,005 and 007 for 3of4.txt",sep="\t")


h_t10=read.table("patterns for  thr_0.1 _subset_3of4_ 2010-11-30 .txt",header=TRUE)
colnames(h_t10)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr10=apply(h_t10[,(1:42)],2,sum)
count_rows_thr10=apply(h_t10[,(1:42)],1,sum)
t10=count_cols_thr10
t10=cbind(t10,names(t10))
colnames(t10)=c("thr010","patterns")
compar_t1t2t5t7=read.table("comparison thr001,002,005 and 007 for 3of4.txt",header=TRUE)
compar_t1t2t5t7t10=merge(compar_t1t2t5t7,t10)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-c(3,5)]
table_comp=merge(compar_t1t2t5t7t10,topo)
table_comp3of4=table_comp[table_comp$isl3_pat_only=="y",]
write.table(table_comp3of4,"comparison thr001,002,005,007 and 010 for 3of4.txt",sep="\t")

#####works until here (30.11.2010)

###########compacting in a single table the results for 3 islands form 3of4 and 3 only (30.11.2010)
table_only3=read.table("comparison thr001,002,005,007 and 010 from 3 islands only.txt", header=TRUE)
table_3of4=read.table("comparison thr001,002,005,007 and 010 for 3of4.txt", header=TRUE)
test1=table_only3$thr001+table_3of4$thr001
table_3tot=table_only3[,c(1,2)]
test1=table_only3$thr001+table_3of4$thr001
table_3tot=cbind(table_3tot,test1)
test2=table_only3$thr002+table_3of4$thr002
table_3tot=cbind(table_3tot,test2)
test5=table_only3$thr005+table_3of4$thr005
table_3tot=cbind(table_3tot,test5)
test7=table_only3$thr007+table_3of4$thr007
table_3tot=cbind(table_3tot,test7)
test10=table_only3$thr010+table_3of4$thr010
table_3tot=cbind(table_3tot,test10)

names(table_3tot)=c("patterns","Nb_Island","thr001","thr002","thr005","thr007","thr010")
write.table(table_3tot,"comparison thr001,002,005,007 and 010 for 3 islands_total.txt", sep="\t")
table_comp3=table_3tot
table_comp4=read.table("comparison thr001,002,005,007 and 010 from 4 islands.txt", header=TRUE)
names(table_comp4)=c("patterns","Nb_Island","thr001","thr002","thr005","thr007","thr010")
table_comp_tot=rbind(table_comp3,table_comp4)
write.table(table_comp_tot,"table_comparison 4 and 3 island_total.txt", sep="\t")

######################################################################################################################################################
###counting patterns for each threshold
sum_tab_subtr=NULL
aux_problm=NULL
aux=NULL
#thres=list()
#compar=NULL

for(y in 1:length(thres))
	{
	matrixCl=thres[[y]]
	all_subtr=list()
	all_nw=vector()
	aux=which(names(matrixCl)==names(thres[y]))
	uniq1=unique(matrixCl[,aux])
	patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
	pat1=row.names(patterns)

	for (x in uniq1)   
		{
		subAux=matrixCl[matrixCl[,aux]==x,] 
		NbIsland=length(unique(subAux[,2])) 
		aux_1=combn(as.vector(subAux[,1]),NbIsland)
		count_tot=0
 		for (j in 1:dim(aux_1)[2]) 
			{
			test=substr(aux_1[,j],9,9)
			if (length(unique(test))==NbIsland)
				{
				count_tot=count_tot+1
				subtr_nw=NULL # restart subtr_nw
				subtr=subset(my3_phylo4,tips.include=aux_1[,j])
				all_subtr[[paste(x,j,sep="_")]]=subtr
				subtr_phylo=as(subtr,"phylo")
				subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
				subtr_phylo$edge.length=NULL
				subtr_newick=write.tree(subtr_phylo)
				subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
				sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
				if (length(sumary_subtr)==2)
					{
					sumary_subtr=c(sumary_subtr, "pb")
					}
				sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
				patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1
			}			
		patterns[pat1==x,pat2=="total"]=count_tot
		} 
	write.table(patterns,paste("patterns for ",names(thres[y]),"_",Sys.Date(),".txt"), sep="\t",row.names=TRUE,col.names=TRUE)
	write.table(sum_tab_subtr,paste("summary of table subtree",names(thres[y]),"_",Sys.Date(),".txt"))
	save(all_subtr,file=paste("list of all subtrees ",names(thres[y]),"_",Sys.Date(),".RData"))
	}
#########works until here!! 22.11.2010

####################################
######not yet tested!!!!############
count_cols_thr=apply(patterns[,(1:42)],2,sum)
count_rows_thr=apply(patterns[,(1:42)],1,sum)
t=count_cols_thr
t=cbind(t,names(t))
colnames(t)=c(names(thres[y]),"patterns")
compar=merge(compar,t)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
table_comp=merge(compar,topo)
#}


###################################
#######version not optimized but working
h_t1=read.table("patterns for thr_001_18112010.txt",header=TRUE)
colnames(h_t1)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr1=apply(h_t1[,(1:42)],2,sum)
count_rows_thr1=apply(h_t1[,(1:42)],1,sum)
t1=count_cols_thr1
t1=cbind(t1,names(t1))
colnames(t1)=c("thr001","patterns")


h_t2=read.table("patterns for thr_002_18112010.txt",header=TRUE)
colnames(h_t2)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr2=apply(h_t2[,(1:42)],2,sum)
count_rows_thr2=apply(h_t2[,(1:42)],1,sum)
t2=count_cols_thr2
t2=cbind(t2,names(t2))
colnames(t2)=c("thr002","patterns")
compar_t1t2=merge(t1,t2)
#compar_t1t2=compar_t1t2[,-2]
#topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
#table_comp=merge(compar_t1t2,topo)

h_t5=read.table("patterns for thr_005_19112010.txt",header=TRUE)
colnames(h_t5)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr5=apply(h_t5[,(1:42)],2,sum)
count_rows_thr5=apply(h_t5[,(1:42)],1,sum)
t5=count_cols_thr5
t5=cbind(t5,names(t5))
colnames(t5)=c("thr005","patterns")
compar_t1t2t5=merge(compar_t1t2,t5)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-5]
table_comp=merge(compar_t1t2t5,topo)
table_comp3=table_comp[table_comp$isl3_pat_only=="y",]
table_comp4=table_comp[table_comp$isl4_pat=="y",]
write.table(table_comp,"comparison thr001,002 and 005 from 4 and 3 islands.txt",sep="\t")
write.table(table_comp4,"comparison thr001,002 and 005 from 4 islands.txt",sep="\t")
write.table(table_comp3,"comparison thr001,002 and 005 from 3 islands only.txt",sep="\t")

h_t7=read.table("patterns for  thr_0.07 _ 2010-11-28 .txt",header=TRUE)
colnames(h_t7)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr7=apply(h_t7[,(1:42)],2,sum)
count_rows_thr7=apply(h_t7[,(1:42)],1,sum)
t7=count_cols_thr7
t7=cbind(t7,names(t7))
colnames(t7)=c("thr007","patterns")
compar_t1t2t5t7=merge(compar_t1t2t5,t7)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-5]
table_comp=merge(compar_t1t2t5t7,topo)
table_comp3=table_comp[table_comp$isl3_pat_only=="y",]
table_comp4=table_comp[table_comp$isl4_pat=="y",]
write.table(table_comp,"comparison thr001,002,005 and 007 from 4 and 3 islands.txt",sep="\t")
write.table(table_comp4,"comparison thr001,002,005 and 007 from 4 islands.txt",sep="\t")
write.table(table_comp3,"comparison thr001,002,005 and 007 from 3 islands only.txt",sep="\t")


h_t10=read.table("patterns for  thr_0.1 _ 2010-11-30 .txt",header=TRUE)
colnames(h_t10)=c("(r,h,(t,m));","(t,(r,(m,h)));","(r,(t,(m,h)));","(t,(m,(h,r)));","(m,(t,(r,h)));","(r,(m,(t,h)));","(m,(r,(t,h)));","(t,(h,(m,r)));","(h,(t,(m,r)));","(r,(h,(t,m)));","(h,(r,(t,m)));","(m,(h,(t,r)));","(h,(m,(t,r)));","(t,(r,m,h));","(r,(m,t,h));","((r,t),(h,m));","(r,t,(m,h));","(r,(h,m));","(t,(m,h));","(m,(r,t,h));","((r,h),(m,t));","((r,h),t,m);","(m,(h,r));","(t,(h,r));","((r,m),(t,h));","(r,m,(t,h));","(m,(t,h));","(r,(t,h));","(h,(t,m,r));","((r,m),t,h);","(h,(m,r));","(t,(m,r));","(h,(t,m));","(r,(t,m));","((r,t),m,h);","(h,(t,r));","(m,(t,r));","(r,m,t,h);","(m,r,h);","(m,h,t);","(h,r,t);","(m,r,t);","total")
count_cols_thr10=apply(h_t10[,(1:42)],2,sum)
count_rows_thr10=apply(h_t10[,(1:42)],1,sum)
t10=count_cols_thr10
t10=cbind(t10,names(t10))
colnames(t10)=c("thr010","patterns")
compar_t1t2t5t7=read.table("comparison thr001,002,005 and 007 from 4 and 3 islands copy.txt",header=TRUE)
compar_t1t2t5t7t10=merge(compar_t1t2t5t7,t10)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-5]
table_comp=merge(compar_t1t2t5t7t10,topo)
table_comp3=table_comp[table_comp$isl3_pat_only=="y",]
table_comp4=table_comp[table_comp$isl4_pat=="y",]
write.table(table_comp,"comparison thr001,002,005,007 and 010 from 4 and 3 islands.txt",sep="\t")
write.table(table_comp4,"comparison thr001,002,005,007 and 010 from 4 islands.txt",sep="\t")
write.table(table_comp3,"comparison thr001,002,005,007 and 010 from 3 islands only.txt",sep="\t")

####changed on the 30.11.2010 with inclusion of thr010

######implementation of frequency events for all other possible combinations of pairs of islands (x,y) with x,y=m,r,h,t
one=NULL
two=NULL
p_tot=NULL
p_total=NULL
p_fav=NULL
pairs=NULL
case4f=NULL
case4t=NULL
case3f=NULL
case3t=NULL
thr_p=NULL

pairs=c("(m,t)","(m,r)","(m,h)","(t,r)","(t,h)","(r,h)")
thr_p=c("thr001","thr002","thr005","thr007","thr010")
table_comp=read.table("table_comparison 4 and 3 island_total.txt",header=TRUE)

for (p in pairs)
	{
	one=substr(p,2,2)
	two=substr(p,4,4)
	inv_p=paste("(",two,",",one,")",sep="")
	p_fav=c(grep(p,table_comp$patterns, fixed=TRUE, value=TRUE),grep(inv_p, table_comp[,1], fixed=TRUE,value=TRUE))
	p_total=grep(one,grep(two,table_comp$patterns, value=TRUE), value=TRUE) 
	for(e in 1:length(p_total))
		{
		hlp=which(table_comp$patterns==p_total[e])
		p_tot=c(p_tot,hlp)
		}
	
	for (n in p_fav)
		{
		if(table_comp[n,]$Nb_Island==4)
			{
			case4f=c(case4f,n)
			}	
		if(table_comp[n,]$Nb_Island==3)
			{
			case3f=c(case3f,n)
			}
		}
	for (n in p_tot)
		{
		if(table_comp[n,]$Nb_Island==4)
			{
			case4t=c(case4t,n)
			}	
		if(table_comp[n,]$Nb_Island==3)
			{
			case3t=c(case3t,n)
			}
		}

	for (t in 1:length(thr_p))
		{
	###for 4 islands
		Nb_fav4_p=sum(as.numeric(as.matrix(table_comp[case4f,(t+2)])))
	print(Nb_fav4_p)
		Nb_tot4_p=sum(as.numeric(as.matrix(table_comp[case4t,(t+2)])))
	assump4_p=matrix(c(Nb_fav4_p,Nb_tot4_p,(Nb_fav4_p/Nb_tot4_p),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases","freq.obs","freq. exp"),p))
		write.table(assump4_p,paste("assump4_",one,two,thr_p[t],"test.txt",sep=""))

	###for 3 islands
		Nb_fav3_p=sum(as.numeric(as.matrix(table_comp[case3f,(t+2)])))
		Nb_tot3_p=sum(as.numeric(as.matrix(table_comp[case3t,(t+2)])))
	assump3_p=matrix(c(Nb_fav3_p,Nb_tot3_p,(Nb_fav3_p/Nb_tot3_p),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),p))
		write.table(assump3_p,paste("assump3_",one,two,thr_p[t],"test.txt",sep=""))
		}
	}
###########works until here (30.11.2010)



#################################################################################################################################### as working before the 30.11.2010
one=NULL
two=NULL
p_tot=NULL
p_total=NULL
p_fav=NULL
pairs=NULL
case4f=NULL
case4t=NULL
case3f=NULL
case3t=NULL
Nb_fav4_p001=NULL
Nb_tot4_p001=NULL
Nb_fav4_p002=NULL
Nb_tot4_p002=NULL
Nb_fav3_p001=NULL
Nb_tot3_p001=NULL
Nb_fav3_p002=NULL
Nb_tot3_p002=NULL
Nb_fav4_p005=NULL
Nb_tot4_p005=NULL
Nb_fav3_p005=NULL
Nb_tot3_p005=NULL
thr_p=NULL


for (p in pairs)
	{
	one=substr(p,2,2)
	two=substr(p,4,4)
	inv_p=paste("(",two,",",one,")",sep="")
	p_fav=c(grep(p,table_comp$patterns, fixed=TRUE),grep(inv_p, table_comp[,1], fixed=TRUE))
	p_total=grep(one,grep(two,table_comp$patterns, value=TRUE), value=TRUE) 
	for(e in 1:length(p_total))
		{
		hlp=which(table_comp$patterns==p_total[e])
		p_tot=c(p_tot,hlp)
		}
	
	for (n in p_fav)
		{
		if(table_comp[n,]$Nb_island==4 & table_comp[n,]$isl4_pat=="y")
			{
			case4f=c(case4f,n)
			}	
		if(table_comp[n,]$Nb_island==3 & table_comp[n,]$isl3_pat_only=="y")
			{
			case3f=c(case3f,n)
			}
		}
	for (n in p_tot)
		{
		if(table_comp[n,]$Nb_island==4 & table_comp[n,]$isl4_pat=="y")
			{
			case4t=c(case4t,n)
			}	
		if(table_comp[n,]$Nb_island==3 & table_comp[n,]$isl3_pat_only=="y")
			{
			case3t=c(case3t,n)
			}
		}

	for (t in 1:length(thr_p))
		{
	###for 4 islands
		Nb_fav4_p=sum(as.numeric(as.matrix(table_comp[case4f,(t+1)])))
	print(Nb_fav4_p)
		Nb_tot4_p=sum(as.numeric(as.matrix(table_comp[case4t,(t+1)])))
	assump4_p=matrix(c(Nb_fav4_p,Nb_tot4_p,(Nb_fav4_p/Nb_tot4_p),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases","freq.obs","freq. exp"),p))
		write.table(assump4_p,paste("assump4_",one,two,t,"test.txt",sep=""))

	###for 3 islands
		Nb_fav3_p=sum(as.numeric(as.matrix(table_comp[case3f,(t+1)])))
		Nb_tot3_p=sum(as.numeric(as.matrix(table_comp[case3t,(t+1)])))
	assump3_p=matrix(c(Nb_fav3_p,Nb_tot3_p,(Nb_fav3_p/Nb_tot3_p),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),p))
		write.table(assump3_p,paste("assump3_",one,two,thr_p[t],"test.txt",sep=""))
		}
	}


########################################################################################################################################################################
#############need to implement for cases (h,(m,t)) on the 19.11.2010
	
### (h,(m,t)) cases
hmt=c(grep("h,(m,t)",table_comp[,1], fixed=TRUE),grep("h,(t,m)", table_comp[,1], fixed=TRUE),grep("(m,t),h",table_comp[,1], fixed=TRUE),grep("(t,m),h",table_comp[,1], fixed=TRUE))
hmt_total=grep("h",grep("t",grep("m",table_comp[,1], value=TRUE), value=TRUE),value=TRUE) 
hmt_tot=NULL
for(d in 1:length(hmt_total))
	{
		hlp=which(table_comp[,1]==hmt_total[d])
		hmt_tot=c(hmt_tot,hlp)
	}

hmt3=NULL
for(h in hmt)
	{
	if (table_comp[h,5]=="3")
		{
		hmt3=c(hmt3,h)
		}
	}
Nb_fav3_hmt001=sum(as.numeric(as.matrix(table_comp[hmt3,3])))





### (r,(h,(t,m))) cases
rhmt=grep("(r,(h,(t,m)))",table_comp[,1], fixed=TRUE)
rhmt_total=table_comp[table_comp[,5]=="y",1]
rmht_tot=NULL
for(c in 1:length(rhmt_total))
	{
		hlp=which(table_comp[,1]==rhmt_total[c])
		rhmt_tot=c(rhmt_tot,hlp)
	}
	



##producing a summary tables for 3 and 4 islands @thr=0.01

assump3_001=matrix(c(Nb_fav3_mt001,Nb_tot3_001,Nb_fav3_hmt001, Nb_tot3_001), nrow=2, ncol=2, dimnames=list(c("favourable cases", "total cases"),c("(m,t)","h,(m,t)")))

assump4_001=matrix(c(Nb_fav4_mt001,Nb_tot4_mt001,Nb_fav4_hmt001, Nb_tot4_hmt001,Nb_fav4_rhmt001,Nb_tot4_rhmt001), nrow=2, ncol=3, dimnames=list(c("favourable cases", "total cases"),c("(m,t)","h,(m,t)","r,(h,(m,t))")))

assump4_002=matrix(c(Nb_fav4_mt002,Nb_tot4_mt002,Nb_fav4_hmt002, Nb_tot4_hmt002,Nb_fav4_rhmt002,Nb_tot4_rhmt002), nrow=2, ncol=3, dimnames=list(c("favourable cases", "total cases"),c("(m,t)","h,(m,t)","r,(h,(m,t))")))



########################################################################################################################


####test if column total and sum of row are equal
#aux_problm2=NULL
#for(m in names(count_rows))
#	{
#		if(patterns[row.names(patterns)==m,43]==count_rows[names(count_rows)==m])
#		print("no problem")
#		else{ print(paste("problem with",m,sep="\t"))
#			aux_problm2=c(aux_problm2,m)}
#	}


#count_cols002=apply(patterns[,(1:42)],2,sum)
#count_rows002=apply(patterns[,(1:42)],1,sum)





##################################################################################################################################################################################################################################################################################################################################################################################################################################################################################################################
##################as implemented to calculate individually thr=0.001,0.02 and 0.05 before the 22.11.2010

island_001=read.table("island3to4_thr_0.01")
island_002=read.table("island3to4_thr_0.02")
island_005=read.table("island3to4_thr_0.05")
island_007=read.table("island3to4_thr_0.07")
island_01=read.table("island3to4_thr_0.1")


##############################
###patterns for threshold=0.01
sum_tab_subtr=NULL

matrixCl=island_001
all_subtr001=list()
uniq1=unique(matrixCl$thr_0.01)
#table_patterns=read.table("table_patterns.txt",header=TRUE,colClasses="character")
patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
pat1=row.names(patterns)
all_nw001=vector()


for (x in uniq1)   
{
	subAux=matrixCl[matrixCl[,6]==x,] 
	NbIsland=length(unique(subAux[,2])) 
	aux_1=combn(as.vector(subAux[,1]),NbIsland)
	count_tot=0
 	for (j in 1:dim(aux_1)[2]) 
	{
		test=substr(aux_1[,j],9,9)
	#print(test)
		if (length(unique(test))==NbIsland)
		{
			count_tot=count_tot+1
			subtr_nw=NULL # restart subtr_nw
			subtr=subset(my3_phylo4,tips.include=aux_1[,j])
			all_subtr001[[paste(x,j,sep="_")]]=subtr
			subtr_phylo=as(subtr,"phylo")
			subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
			subtr_phylo$edge.length=NULL
			subtr_newick=write.tree(subtr_phylo)
			subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
			sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
			if (length(sumary_subtr)==2)
			{
				sumary_subtr=c(sumary_subtr, "pb")
			}
			#print(sumary_subtr)
			sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
			#all_nw001=c(all_nw001,subtr_nw)
			patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1
		}			
	} 
	#print(paste("Cluster",x,"done\n",sep=" "))
	patterns[pat1==x,pat2=="total"]=count_tot
}
write.table(patterns,"patterns for thr_001_18112010.txt", sep="\t",row.names=TRUE,col.names=TRUE)
write.table(sum_tab_subtr,"summary of table subtree.txt")
save(all_subtr001,file="list of all subtrees thr_001_18112010.RData")

####counting total events for each pattern (columns) and for each cluster (rows)
count_cols_thr1=apply(patterns[,(1:42)],2,sum)
count_rows_thr1=apply(patterns[,(1:42)],1,sum)
t1=count_cols_thr1
t1=cbin(t1,names(t1))
colnames(t1)=c("thr001","patterns")

####test if column total and sum of row are equal
aux_problm=NULL
for(m in names(count_rows))
	{
	if(patterns[row.names(patterns)==m,43]==count_rows[names(count_rows)==m])
	print("no problem")
		else{ print(paste("problem with",m,sep="\t"))
			aux_problm=c(aux_problm,m)}
	}
	
	
##############################	
#####applied to thr=0.02
sumary_subtr=NULL
sum_tab_subtr=NULL

matrixCl=island_002
all_subtr002=list()
uniq1=unique(matrixCl$thr_0.02)
patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
pat1=row.names(patterns)

for (x in uniq1)   
{
	subAux=matrixCl[matrixCl[,6]==x,] 
	NbIsland=length(unique(subAux[,2])) 
	aux_1=combn(as.vector(subAux[,1]),NbIsland)
	count_tot=0
 	for (j in 1:dim(aux_1)[2]) 
	{
		test=substr(aux_1[,j],9,9)
	#print(test)
		if (length(unique(test))==NbIsland)
		{
			count_tot=count_tot+1
			subtr_nw=NULL # restart subtr_nw
			subtr=subset(my3_phylo4,tips.include=aux_1[,j])
			all_subtr002[[paste(x,j,sep="_")]]=subtr
			subtr_phylo=as(subtr,"phylo")
			subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
			subtr_phylo$edge.length=NULL
			subtr_newick=write.tree(subtr_phylo)
			subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
			sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
			if (length(sumary_subtr)==2)
			{
				sumary_subtr=c(sumary_subtr, "pb")
			}
			#print(sumary_subtr)
			sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
			#all_nw001=c(all_nw001,subtr_nw)
			patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1
		}			
	} 
	#print(paste("Cluster",x,"done\n",sep=" "))
	patterns[pat1==x,pat2=="total"]=count_tot
}
write.table(patterns,"patterns for thr_002_18112010.txt", sep="\t",row.names=TRUE,col.names=TRUE)
write.table(sum_tab_subtr,"summary of table subtree002.txt")
save(all_subtr002,file="list of all subtrees thr_002_18112010.RData")


count_cols_thr2=apply(patterns[,(1:42)],2,sum)
count_rows_thr2=apply(patterns[,(1:42)],1,sum)

t2=count_cols_thr2
t2=cbin(t2,names(t2))
colnames(t2)=c("thr002","patterns")
compar_t1t2=merge(t1,t2)
compar_t1t2=compar_t1t2[,-2]
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
table_comp=merge(compar_t1t2,topo)


##############################
#####applied to thr=0.05
sumary_subtr=NULL
sum_tab_subtr=NULL
subAux=NULL
NbIsland=NULL
aux_1=NULL
test=NULL
subtr=NULL
subtr_phylo=NULL
subtr_newick=NULL
subtr_nw=NULL
patterns=NULL
matrixCl=NULL
t5=NULL
compar_t1t2t5=NULL

matrixCl=island_005
all_subtr005=list()
uniq1=unique(matrixCl$thr_0.05)
patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
pat1=row.names(patterns)

for (x in uniq1)   
	{
	subAux=matrixCl[matrixCl$thr_0.05==x,] 
	NbIsland=length(unique(subAux[,2])) 
	aux_1=combn(as.vector(subAux[,1]),NbIsland)
	count_tot=0
 	for (j in 1:dim(aux_1)[2]) 
		{
		test=substr(aux_1[,j],9,9)
		if (length(unique(test))==NbIsland)
			{
			count_tot=count_tot+1
			subtr_nw=NULL # restart subtr_nw
			subtr=subset(my3_phylo4,tips.include=aux_1[,j])
			all_subtr005[[paste(x,j,sep="_")]]=subtr
			subtr_phylo=as(subtr,"phylo")
			subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
			subtr_phylo$edge.length=NULL
			subtr_newick=write.tree(subtr_phylo)
			subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
			sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
			if (length(sumary_subtr)==2)
				{
				sumary_subtr=c(sumary_subtr, "pb")
				}
			sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
			patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1
			}			
		} 
		patterns[pat1==x,pat2=="total"]=count_tot
	}
write.table(patterns,"patterns for thr_005_19112010.txt", sep="\t",row.names=TRUE,col.names=TRUE)
write.table(sum_tab_subtr,"summary of table subtree005.txt")
save(all_subtr005,file="list of all subtrees thr_005_19112010.RData")

count_cols_thr5=apply(patterns[,(1:42)],2,sum)
count_rows_thr5=apply(patterns[,(1:42)],1,sum)


t5=count_cols_thr5
t5=cbind(t5,names(t5))
colnames(t5)=c("thr005","patterns")
compar_t1t2t5=merge(compar_t1t2,t5)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-5]
table_comp=merge(compar_t1t2t5,topo)

#########################################################################################################################################################################################################################################
####code as it worked on 19.11.2010 when table_comp had only thr=0.01 and 0.02#####
##############################################################################code with access to columns made with $ and tested to create mt005#######
##### gives row numbers of favourable cases for our hypothesis
### (m,t) cases
mt=NULL
mt_total=NULL
mt_tot=NULL
case4f=NULL
case3f=NULL
case4t=NULL
case3t=NULL
Nb_fav4_mt001=NULL
Nb_tot4_mt001=NULL
Nb_fav4_mt002=NULL
Nb_tot4_mt002=NULL
Nb_fav3_mt001=NULL
Nb_tot3_mt001=NULL
Nb_fav3_mt002=NULL
Nb_tot3_mt002=NULL
Nb_fav4_mt005=NULL
Nb_tot4_mt005=NULL
Nb_fav3_mt005=NULL
Nb_tot3_mt005=NULL


mt=c(grep("(m,t)",table_comp$patterns, fixed=TRUE),grep("(t,m)", table_comp[,1], fixed=TRUE))
mt_total=grep("t",grep("m",table_comp$patterns, value=TRUE), value=TRUE) 
for(e in 1:length(mt_total))
	{
	hlp=which(table_comp$patterns==mt_total[e])
	mt_tot=c(mt_tot,hlp)
	}

#####separates favourable patterns into 3 or 4 islands ones that can be taken as possible events
for (n in mt)
	{
	if(table_comp[n,]$Nb_island==4 & table_comp[n,]$isl4_pat=="y")
		{
		case4f=c(case4f,n)
		}	
	if(table_comp[n,]$Nb_island==3 & table_comp[n,]$isl3_pat_only=="y")
		{
		case3f=c(case3f,n)
		}
	}
	
#####separates possible patterns into 3 or 4 islands ones that can be taken as possible events
for (n in mt_tot)
	{
	if(table_comp[n,]$Nb_island==4 & table_comp[n,]$isl4_pat=="y")
		{
		case4t=c(case4t,n)
		}	
	if(table_comp[n,]$Nb_island==3 & table_comp[n,]$isl3_pat_only=="y")
		{
		case3t=c(case3t,n)
		}
	}

###for 4 islands
##thr=0.01
Nb_fav4_mt001=sum(as.numeric(as.matrix(table_comp[case4f,2])))
Nb_tot4_mt001=sum(as.numeric(as.matrix(table_comp[case4t,2]
)))
assump4_mt001=matrix(c(Nb_fav4_mt001,Nb_tot4_mt001,(Nb_fav4_mt001/Nb_tot4_mt001),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases","freq.obs","freq. exp"),"(m,t)"))
#write.table(assump4_mt001,"assump4_mt001.txt")

##thr=0.02
Nb_fav4_mt002=sum(as.numeric(as.matrix(table_comp[case4f,3])))
Nb_tot4_mt002=sum(as.numeric(as.matrix(table_comp[case4t,3])))
assump4_mt002=matrix(c(Nb_fav4_mt002,Nb_tot4_mt002,(Nb_fav4_mt002/Nb_tot4_mt002),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
#write.table(assump4_mt002,"assump4_mt002.txt")

##thr=0.05
Nb_fav4_mt005=sum(as.numeric(as.matrix(table_comp[case4f,4])))
Nb_tot4_mt005=sum(as.numeric(as.matrix(table_comp[case4t,4])))
assump4_mt005=matrix(c(Nb_fav4_mt005,Nb_tot4_mt005,(Nb_fav4_mt005/Nb_tot4_mt005),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases","freq.obs","freq. exp"),"(m,t)"))
write.table(assump4_mt005,"assump4_mt005.txt")



###for 3 islands
##thr=0.01
Nb_fav3_mt001=sum(as.numeric(as.matrix(table_comp[case3f,2])))
Nb_tot3_mt001=sum(as.numeric(as.matrix(table_comp[case3t,2])))
assump3_mt001=matrix(c(Nb_fav3_mt001,Nb_tot3_mt001,(Nb_fav3_mt001/Nb_tot3_mt001),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
#write.table(assump3_mt001,"assump3_mt001.txt")

##thr=0.02
Nb_fav3_mt002=sum(as.numeric(as.matrix(table_comp[case3f,3])))
Nb_tot3_mt002=sum(as.numeric(as.matrix(table_comp[case3t,3])))
assump3_mt002=matrix(c(Nb_fav3_mt002,Nb_tot3_mt002,(Nb_fav3_mt002/Nb_tot3_mt002),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
#write.table(assump3_mt002,"assump3_mt002.txt")

##thr=0.05
Nb_fav3_mt005=sum(as.numeric(as.matrix(table_comp[case3f,4])))
Nb_tot3_mt005=sum(as.numeric(as.matrix(table_comp[case3t,4])))
assump3_mt005=matrix(c(Nb_fav3_mt005,Nb_tot3_mt005,(Nb_fav3_mt005/Nb_tot3_mt005),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
write.table(assump3_mt005,"assump3_mt005.txt")


##### gives row numbers of favourable cases for our hypothesis
### (m,t) cases
mt=NULL
mt_total=NULL
mt_tot=NULL
case4f=NULL
case3f=NULL
case4t=NULL
case3t=NULL
Nb_fav4_mt001=NULL
Nb_tot4_mt001=NULL
Nb_fav4_mt002=NULL
Nb_tot4_mt002=NULL
Nb_fav3_mt001=NULL
Nb_tot3_mt001=NULL
Nb_fav3_mt002=NULL
Nb_tot3_mt002=NULL


mt=c(grep("(m,t)",table_comp[,1], fixed=TRUE),grep("(t,m)", table_comp[,1], fixed=TRUE))
mt_total=grep("t",grep("m",table_comp[,1], value=TRUE), value=TRUE) # all cases where t and m appear
mt_tot=NULL
for(e in 1:length(mt_total))
	{
		hlp=which(table_comp[,1]==mt_total[e])
		mt_tot=c(mt_tot,hlp)
	}
	
	
#####separates favourable patterns into 3 or 4 islands ones that can be taken as possible events
for (n in mt)
	{
	if(table_comp[n,4]==4 & table_comp[n,5]=="y")
		{
			case4f=c(case4f,n)
		}	
	if(table_comp[n,4]==3 & table_comp[n,6]=="y")
		{
			case3f=c(case3f,n)
		}
	}
	
#####separates possible patterns into 3 or 4 islands ones that can be taken as possible events
for (n in mt_tot)
	{
	if(table_comp[n,4]==4 & table_comp[n,5]=="y")
		{
			case4t=c(case4t,n)
		}	
	if(table_comp[n,4]==3 & table_comp[n,6]=="y")
		{
			case3t=c(case3t,n)
		}
	}

###for 4 islands
##thr=0.01
Nb_fav4_mt001=sum(as.numeric(as.matrix(table_comp[case4f,2])))
Nb_tot4_mt001=sum(as.numeric(as.matrix(table_comp[case4t,2]
)))
assump4_mt001=matrix(c(Nb_fav4_mt001,Nb_tot4_mt001,(Nb_fav4_mt001/Nb_tot4_mt001),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases","freq.obs","freq. exp"),"(m,t)"))
write.table(assump4_mt001,"assump4_mt001.txt")

##thr=0.02
Nb_fav4_mt002=sum(as.numeric(as.matrix(table_comp[case4f,3])))
Nb_tot4_mt002=sum(as.numeric(as.matrix(table_comp[case4t,3]
)))
assump4_mt002=matrix(c(Nb_fav4_mt002,Nb_tot4_mt002,(Nb_fav4_mt002/Nb_tot4_mt002),0.20), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
write.table(assump4_mt002,"assump4_mt002.txt")


###for 3 islands
##thr=0.01
Nb_fav3_mt001=sum(as.numeric(as.matrix(table_comp[case3f,2])))
Nb_tot3_mt001=sum(as.numeric(as.matrix(table_comp[case3t,2]
)))
assump3_mt001=matrix(c(Nb_fav3_mt001,Nb_tot3_mt001,(Nb_fav3_mt001/Nb_tot3_mt001),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
write.table(assump3_mt001,"assump3_mt001.txt")


##thr=0.02
Nb_fav3_mt002=sum(as.numeric(as.matrix(table_comp[case3f,3])))
Nb_tot3_mt002=sum(as.numeric(as.matrix(table_comp[case3t,3])))

assump3_mt002=matrix(c(Nb_fav3_mt002,Nb_tot3_mt002,(Nb_fav3_mt002/Nb_tot3_mt002),0.33), nrow=4, ncol=1, dimnames=list(c("favourable cases", "total cases", "freq.obs","freq. exp"),"(m,t)"))
write.table(assump3_mt002,"assump3_mt002.txt")

